Towards active microfluidics: Interface turbulence in thin liquid films with floating 

molecular machines 
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Thin liquid films with floating active protein machines are considered. Cyclic mechanical motions 
within the machines, representing microscopic swimmers, lead to molecular propulsion forces applied 
to the air-liquid interface. We show that, when the rate of energy supply to the machines exceeds a 
threshold, the flat interface becomes linearly unstable. As the result of this instability, the regime of 
interface turbulence, characterized by irregular traveling waves and propagating machine clusters, 
is established. Numerical investigations of this nonlinear regime are performed. Conditions for the 
experimental observation of the instability are discussed. 



I. INTRODUCTION 

Molecular machines are protein molecules which can 
transform chemical energy into ordered internal mechan- 
ical motions. The classical examples of such machines 
are molecular motors kinesin and myosin, where internal 
mechanical motions are used to transport cargo along 
microtubules and filaments. Many enzymes operate as 
machines, using internal conformational motions to facil- 
itate chemical reactions. Other kinds of machines, oper- 
ating as ion pumps or involved in genetic processes, are 
also known. Moreover, artificial nonequilibrium nanodc- 
vices, similar to protein machines, are being developed 

a. 

The cycle of a protein machine is typically powered by 
chemical energy supplied with an ATP molecule. When 
an ATP molecule binds to a protein machine, the initial 
protein- ATP complex is out of equilibrium and the pro- 
cess of ordered conformational relaxation of this complex 
towards its equilibrium state starts. When this state is 
reached, ATP is converted into ADP and this product 
molecule leaves the complex. The emerging free pro- 
tein molecule is out of equilibrium and another confor- 
mational relaxation process, returning the free protein 
to its original equilibrium state begins. Thus, the cycle 
consists of two conformational motions. It is important 
to note that the forward and back conformational mo- 
tions inside each cycle do not coincide - they correspond 
to relaxation processes in two different physical objects, 
i.e. the free protein and the protein- ATP complex Q. 

Recently, much attention has been attracted to mi- 
croscale swimmers operating at low Reynolds numbers. 
Generally, it can be shown that any physical object, cycli- 
cally changing its shape in such a way that the inter- 
nal forward and back motions are different, propels itself 
through the liquid [!, Q|. Elementary models of such 
swimmers, constructed by joining together mobile links 
H, IE] or by connecting a few spheres by several actively 
deformable links @, 0, [1], have been considered. Typ- 
ically, the concept of molecular swimmers is applied to 
explain active motion of bacteria and other microorgan- 
isms. However, individual macromolecules that operate 



actively as machines are also capable of self-motion [l|, |j| 
If a swimmer is attached to some support, preventing 
its translational motion, it exhibits force acting on the 
mechanical support. The generated force is equal to the 
stall force which needs to be applied to a swimmer to pre- 
vent its translational motion. When many swimmers are 
attached to a distributed support, their collective opera- 
tion produces mechanical pressure acting on the support. 
This situation is, for example, encountered for molecular 
machines representing active protein inclusions (such as 
ion pumps) in biological membranes. It has been shown 
that the combined effects of the mechanical forces gen- 
erated by the machines and their lateral motions inside 
the membrane can lead to membrane instabilities [lj| . 

The aim of this paper is to investigate a different situ- 
ation involving molecular machine populations. We con- 
sider hypothetical machines which arc floating on top of 
a thin liquid layer and thus represent active surfactants 
(FigOJ. Because the considered surfactant molecules are 
microswimmers, they generate local pressure applied to 
the liquid-air interface and proportional to their local 
concentration. Additionally, they are subject to surface 
diffusion and, as surfactants, also modify local surface 
tension. Our main analytical result, supported by numer- 
ical simulations, is that, if the surface density of machines 
is sufficiently high and if the rate of supply of chemical 
energy to the machine population exceeds a threshold, 
the equilibrium flat interface becomes unstable and hy- 
drodynamical flows inside the liquid layer spontaneously 
develop. This instability is accompanied by spatial re- 
distribution of floating machines and their collective mo- 
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FIG. 1: (Color online) Sketch of a thin film with floating 
molecular machines generating propulsion forces. 
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tions over the surface. It leads to the emergence of a 
special kind of surface turbulence. Such phenomena can 
provide the basis for active microfluidics, where hydro- 
dynamical motions in liquid layers are induced and con- 
trolled by molecular machines located at the surface. 

The next two sections are devoted to the derivation 
and the stability analysis of a model of such molecular 
swimmers attached to the surface of a thin liquid film. A 
collection of numerical results is presented in section IV. 
Finally, a model and numerical results are discussed, and 
some estimates of the forces of such molecules are given. 
The possible experimental realization is also discussed. 



II. FORMULATION OF THE MODEL 

Population of identical molecular machines floating on 
the surface of a thin liquid film will be considered. These 
machines perform cycles of conformational changes (with 
a characteristic time t c ) enabled by the supply of energy 
through ATP molecules present in the liquid. We are 
interested in macroscopic collective effects and do not 
specify the details of machine operation. It will be as- 
sumed that, on the average, each machine generates the 
force which is applied to the air-liquid interface. This 
average force acts along the normal direction (the forces 
in the lateral direction vanish after averaging because of 
the rotational diffusion). Each machine cycle is initiated 
by binding of an ATP molecule and the average force is 
proportional to the cycle frequency. Assuming that bind- 
ing of ATP molecules follows the Michaelis-Menten law, 
the average force per machine molecule is 



a) 



/ = fo 



[ATP] 



Katp + [ATP] ■ 



(1) 



Here, fo is the maximal force under the ATP saturation 
conditions, [ATP] denotes the ATP concentration in the 
liquid, and Katp is the characteristic concentration at 
which saturation begins. 

Because of the cycles, the floating machines produce 
additional pressure acting on the air-liquid interface. 
This pressure p m is proportional to the local surface con- 
centration of the machines and the average force gener- 
ated by an individual machine, i.e. p m = f m c, where 
f m . is the force per Mol of molecules (f m = JNa where 
Na is the Avogadro number). Note that the machines 
are still sufficiently well separated one from another and 
possible synchronization effects of their cycles (cf. [ll|) 
are therefore neglected. 

Before proceeding to the detailed formulation of the 
model, we want to outline the origin of the expected sur- 
face instability. Suppose that the average force, gener- 
ated by machines is directed upwards and therefore the 
machines are pulling the liquid up in the vertical direc- 
tion (Fig. 2A). If machine concentration is increased in 
some region, the pulling pressure is higher in this re- 
gion, inducing local rise in the liquid film thickness. This 
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FIG. 2: (Color online) Mechanisms of (a) instability develop- 
ment for upwards directed propulsion forces, and of (b) insta- 
bility damping for downwards directed propulsion forces. 



however leads to lateral hydrodynamical flows which are 
directed inwards and bring even more machines into the 
region. As a result, the positive feedback, responsible for 
the instability, is established. Note that surface diffusion 
of floating machines and capillary forces are acting in the 
opposite direction, suppressing the instability. Thus, it 
is observed only if the average force / generated by a ma- 
chine is large enough, setting a threshold in terms of the 
energy supply rate (i.e., of the ATP concentration in the 
solution). The instability is not possible when the force 
is downwards directed, thus inducing local depressions of 
the liquid layer (Fig. 2B), Then, in contrast to the pre- 
vious case, hydrodynamical flows remove machines from 
the depression, restoring the equilibrium flat film. 

Temporal evolution of the surface concentration field c 
is described by the equation : 



d t c + V(uc) = o?V 2 c, 



(2) 



where u is the lateral fluid velocity and d is the surface 
diffusion. 

Floating proteins represent surfactants and reduce the 
surface tension of the interfaces (see [l2j )■ To take this 
effect into account, we assume that the surface tension 
coefficient a decreases linearly with the machine concen- 
tration, 



Co 



a c c. 



(3) 



Hydrodynamical flows, induced by the gradients in ma- 
chine concentration, should be further considered. We 
assume that the liquid layer is so thin that the lubrica- 
tion approximation, typically employed in microfluidics 
[lH, is justified. As shown in Appendix A, the evolution 
equation for the local hight h of the interface has the 
form 
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(4) 



where the local pressure is p = — erV 2 /i — p m and \x is the 
viscosity of the fluid. Determining the lateral flow veloc- 
ity at the interface (see Appendix A) and substituting 
it into the evolution equation for the surface concentra- 
tion, a closed set of two partial differential equations is 
obtained. 

Explicitly, the considered dynamics of thin films with 
active surfactants is described by equations: 



d t h = V 



r/j3 
3 
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V (cr V 2 /i - a c cV 2 h + f m c) 
Vc 



(5) 



rate of energy supply to the system. Finally, D is the 
dimcnsionless diffusion coefficient of floating machines. 

Comparing the last terms in equation ([8]), we notice 
that spreading of floating machines, induced by changes 
in the surface tension, has the same functional form as 
surface diffusion. If relative variations of the film thick- 
ness and machine concentration are small (H ~ C ~ 1), 
the effective diffusion coefficient of this process is A. Our 
estimates below in Sect. V indicate that the genuine dif- 
fusion constant D of floating machines is typically much 
smaller than the effective diffusion constant A. Having 
this in mind, we retain the terms including the coefficient 
D in our analytical investigations, but put D = when 
numerical simulations of the model are performed. Re- 
sults with non-zero diffusion coefficient have been also ob- 
tained but they are very close to the results with D = 0. 

As shown in Appendix A, the dimcnsionless horizontal 
(U) and vertical (W) velocities of the liquid film at hight 
Z = z/ho and horizontal spatial location X can be found 



d t c = — V 
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U(X, Z) = -AVCZ 



(9) 



f / 



( BVC + V ((1 - AC) V 2 H) ) (2HZ - Z 2 ), 



Note that we have assumed that the active surfactant 
is insoluble [HE!]. 

For subsequent analysis, it is convenient to write these 
equations in the dimcnsionless form. The (equilibrium) 
liquid layer thickness ho will be used as the length unit, 
time will be measured in units of /x/io/co, and local con- 
centration c in units of the equilibrium machine concen- 
tration co- Changing the variables as T = (cro/fiho)t, 
H = h/ho, X = x/ho and C = c/cq, we obtain 



d T H = -iv (ff 3 V [(I - AC) (V 2 H) + BC] 
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(7) 



d T C = -lv 



CH 2 \I ((I - AC) (V 2 iJ) + BC) 



AV 



CHVC 



DV A C. 



(8) 



The new model equations include only three dimen- 
sionlcss parameters 



W(X, Z) = --AV 2 CZ 2 (10) 
+ i (HZ 2 - ^) V 2 (BC + (I - AC) V 2 H) 



Z 2 VH ( (B - AV 2 H)VC + (I - AC)W 2 H 



when the fields C and H are known. Both velocities 
vanish at the bottom of the liquid film, at Z — 0, because 
of the no-penetration and no-slip boundary conditions 
imposed there. 



III. LINEAR STABILITY ANALYSIS 

The model always has the stationary uniform state 
C = H = 1. To perform the linear stability anal- 
ysis of this state, small perturbations H = 1 + AiJ 
exp(iKX + ST) and C = 1 + AC exp(iKX + ST) in 
the form of plane waves with the wavenumber K are in- 
troduced. Linearizing evolution equations with respect 
to small perturbations AH and AC and solving the lin- 
earized equations, growth rates S = S(K) of such per- 
turbations are obtained. 
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The parameter A specifies the characteristic strength of 
floating machines as the surfactant species (decreasing 
the local surface tension of the interface) . The parameter 
B specifies the magnitude of the pressure generated by 
the cycling molecular machines; it is controlled by the 



S(K) = -1( A -^+d\ K 2 - 1(1 - A)K± (11) 
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FIG. 3: (Color online) Real (solid line) and imaginary (light 
line) parts of the growth rate S as functions of the wavenum- 
ber K of the perturbation for (a) B = 0.38, (b) B = 0.4 and 
(c) B = 0.42. Other parameters are A = 0.2 and D = 0. 



It can be easily checked that, in absence of the energy 
supply (B = 0), both rates are real and negative, so that 
the equilibrium flat film is stable as should be expected. 
When the parameter B is increased, the instability de- 
velops at B = B c where 



B c = 2{A + D). 



(12) 



Above the instability threshold, traveling plane waves 
with the wavenumbers K near K = K C (B), 



K C (B) 
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2 1-A 
and the frequencies near f2 = Q C (B), 



(13) 



Q C (B) « ^V2B—3A {B i B j 



3/2 



(14) 



are growing. The fastest growing mode with K = K c and 
Q = Q c is characterized by the growth rate 



Re[5 c (B)] 



3 (B - B c ) 



I- A 
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(15) 



Figure 3 shows dependences Re [S(K )] and Im )] 
at three values of the parameter B below and above the 
instability boundary. While the fastest growing mode 
at B > B c is always oscillatory, with Im [S'(-K')] 7^ 0, 
above the instability boundary the system also always 
has some standing growing modes with the wavenum- 
bers K smaller than K c . As the boundary is approached 
from above, both the wavenumber K c and the frequency 
fi c decrease and vanish at B = B c . The region with two 
real modes S\(K) and S2{K) always lies below K c ; it 
shrinks and vanishes at B = B c . Similar long- wavelength 
instabilities have previously been discussed for other con- 
servative systems (see (l5|). 



IV. NUMERICAL INVESTIGATIONS OF THE 
NONLINEAR REGIME 

To investigate the behavior of the system in the non- 
linear regime above the instability onset, numerical sim- 
ulations have been performed. Equations (|7l8p were in- 
tegrated using the semi-implicit method (see appendix 
B) for a one-dimensional system using periodic bound- 
ary conditions. As the initial condition, the flat inter- 
face with a uniform machine distribution was chosen and 
small random initial perturbations were applied. 

Our main observation is that the instability develop- 
ment results in the emergence of a complex spatiotcm- 
poral regime which can be described as interface turbu- 
lence. This regime is characterized by spontaneous ap- 
pearance of traveling machine clusters (i.e., of spatial re- 
gions where the local machine concentration is increased) 
and of the accompanying local interface modulations. 

Figure(j4]) gives an illustration of the turbulent regimes 
observed at different different distances from the insta- 
bility threshold. The local film thickness H is displayed 
here in gray scale depending on the spatial coordinate 
(the vertical axis) and time. In the slightly supercritical 
regime at B = 0.42, the transient development of stand- 
ing waves is first observed, which is then followed by the 
emergence of an irregular pattern of traveling and collid- 
ing waves. At larger deviations from the critical point, 
the transients arc faster and the irregular wave dynamics 
appears soon after the instability onset. The character- 
istic spatial scale of the turbulence decreases with the 
control parameter B, consistent with the predictions of 
the linear stability analysis. The characteristic velocity 
of traveling waves is also growing with B. 

Figure [5] shows snapshots of computed turbulent pat- 
terns at different deviations from the instability thresh- 
old. Both the interface profiles and the corresponding 
concentration distributions are presented here. Again, a 
decrease of the characteristic wavelength of the irregular 
spatiotemporal patterns under an increase of the control 
parameter B can be noticed. Moreover, we see that the 
characteristic amplitude of the waves grows with B. To 
illustrate the directions of wave propagation, arrows are 
placed in the insets of this figure, where examples of col- 
liding and traveling waves are shown. These arrows indi- 
cate the directions in which the respective concentration 
and interface profile maxima are shifting at the next time 
moment. Note that the motions of machine clusters (i.e, 
of the local concentration maxima) are typically guiding 
the motions of surface bumps. 

Temporal transients leading to turbulent patterns arc 
characterized in Fig. [6l To construct it, we have de- 
termined the maximum and the minimum values H max 
and flmin of the film thickness as function of time, start- 
ing from the initial moment. Their difference AH(T) = 
H max (T) — H min (T) (where T is the dimensionless time) 
can be chosen to describe the amplitude of the develop- 
ing patterns. As seen in the two top panels in Fig. [51 this 
amplitude first grows (exponentially) and the undergoes 
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FIG. 4: (Color online) Interface turbulence. Space-time diagrams display the local film thickness H (thicker regions correspond 
to bright and thinner regions to dark colors) depending on the spatial coordinate and time for : (a) B — 0.42, (b) B = 0.44, (c) 
B = 0.46 and (d) B = 0.48. Other system parameters are A — 0.2 and D — 0. Numerical integrations for a one-dimensional 
system of length Lq = 512 (with Aa; = 1) and the total time interval of To = 5 ■ 10 5 (with At = 0.01). 



saturation. The transient is shorter for the larger devi- 
ation from the threshold (B = 0.46) and the final mean 
amplitude of the turbulent pattern is also then larger. 
To estimate the characteristic wavenumber Kq(T) of 
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FIG. 5: (Color online) Snapshots of the concentration dis- 
tributions (light lines) and of the respective interface profiles 
(black lines) for (a) B = 0.43, (b) B = 0.46 and (c) B = 0.49. 
Other parameters are A — 0.2 and D = 0. Insets on the bot- 
tom show the collision of two waves (d) and the motion of a 
single wave (e). Arrows in the insets indicate the directions 
of motion of the distribution maxima. 



the developing spatial patterns, their Fourier transforms 
were computed and the positions of dominant maxima 
in the spatial power spectra were determined at different 
time moments. As seen in the two bottom panels in 
Fig. [6j these wavenumbers do not significantly change 
during the transients. The characteristic wavelengths of 
the developed turbulent patterns are therefore not much 
different from those of the critical modes. 

Figure [7] shows temporal dependences of the bight 
H(T) and the concentration C(T) at a fixed point of 
the system in the final turbulent state. The characteris- 
tic time of the oscillations is clearly different for the two 
values of the control parameter. Both properties fluctu- 
ate around some mean values. The fluctuations are more 
rapid farther away from the instability threshold. 

Finally, Fig. [H] displays the dependence of several se- 
lected statistical properties of the final turbulent state 
on the deviation AB = B — B c from the critical point 
B c . The characteristic wavenumbers Kq and amplitudes 
AH have been computed as described above. Time aver- 
ages of these properties and statistical dispersion of the 
data are shown in Figs. [8] The characteristic transient 
times T c (Fig. [Sfc) have been estimated by fitting the 
initial computed time dependence for AH{T) to the ex- 
ponential law, AH(T) ~ exp(T/T c ). The characteristic 
frequency fio of the patterns (Fig. [SJi) is estimated by 
computing their temporal power spectra and determining 
the positions of the dominant maxima. 

Solid curves in Fig. [8] give fits of the simulation data 
to the power laws predicted by the linear stability anal- 
ysis. According to equations (JTHJ) and 1(13]). K c — AB 1 / 2 
and S1 c ~ AB 3 ^ 2 . The transient time is determined by 
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FIG. 6: Initial evolution of the pattern amplitudes AH(T) 
and of the characteristic wavenumbers Kq of the developing 
patterns for two different values B = 0.43 and B — 0.46 of 
the control parameter. Other parameters are A = 0.2 and 
D = 0. The system size is Lq = 512. 



the rate of growth as T c = 1/Rc [S C (B)] and, according 
to equation (fT5)) . we expect that T~ l ~ AB 2 . For the 
mean amplitude AH of the patterns, the quadratic fit 
AH ~ AB 1 / 2 has been applied. We see that the sta- 
tistical properties of nonlinear patterns in the developed 
turbulence regime are still in good agreement with the 
respective predictions based on the linear stability anal- 
ysis. 



V. DISCUSSION AND CONCLUSIONS 

Can predicted interface instabilities be experimentally 
observed? This question has both chemical and physical 
aspects. On the chemical side, it is known that proteins 
may indeed represent surfactants and thus float at the 
air- water interface (see. e.g., [H|). Although we cannot 
give here a specific example, it seems plausible that some 
protein machines also belong to this class. Moreover, 
other protein machines, including molecular motors, can 
probably be made floating by chemical modification, i.e. 
by attaching to them a hydrophobic group. 

The physical question is whether the propulsion forces 
generated by individual protein machines would be suf- 
ficient to induce the considered interface instability. Ac- 
cording to equation (Tl12|) , the instability is reached when 
B > B c with B c = 2(A + D). Taking into account the 
definitions of dimcnsionlcss properties A, B and D, the 
instability condition implies that the propulsion force /, 
generated by a single machine, must exceed the threshold 



FIG. 7: Time dependences of the film hight H(T) and concen- 
tration C(T) at a fixed spatial point for two different values of 
the control parameter B. Other parameters and the system 
length are the same as in Fig|4] 



where ho is the film thickness, Co is the surface concen- 
tration of floating machines, d is their surface diffusion 
constant, \x is the liquid viscosity, Na is the Avogadro 
number, and er c specifics the surfactant capacity of the 
considered biomolcculcs (the surface tension coefficient a 
depends as a = <jq — a c c on their concentration c). 

For numerical order-of-magnitude estimates, we choose 
fi ps 1CP 2 g cm _1 s _1 (the viscosity of water), d ss 
10 cm 2 s _1 (characteristic diffusion constant of large 
biomolecules in water solutions |l6IO. Based on the ex- 
perimental data given in Ref. [l2|, we take furthermore 
cr c « 10 7 g cm 2 s~ 2 mMol" x . 

Both terms in equation (|16p arc inversely proportional 
to the liquid film thickness ho and the smaller critical 
forces are therefore expected for the thicker layers. Note, 
however, that this result is based on the lubrication ap- 
proximation for thin liquid films, requiring that the film 
thickness is much smaller than the characteristic wave- 
length of the flow patterns. As we have seen, the charac- 
teristic wavelength of the interface turbulence in the con- 
sidered system depends on the distance from the critical 
point (cf. equation Q13p ) and, in principle, can be arbi- 
trarily large sufficiently close to the instability threshold. 
Taking into account experimental limitations making it 
very difficult to work too close to the threshold, we choose 
however ho = 1mm as the characteristic maximum film 
thickness. 

Concentration Co of floating proteins enters only into 
the second term in equation (|16p . This term, depending 
on surface diffusion of floating machines, can be neglected 
in comparison with the first term, if the condition Co ^> Cd 
with 



U ~ N A \h c Q h 2 



(16) 



a c ho 



(17) 
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FIG. 8: Statistical properties of the interface turbulence at 
different deviations from the instability threshold: (a) char- 
acteristic wavenumber, (b) characteristic frequency, (c) char- 
acteristic inverse transient time and (d) mean amplitude of 
the patterns. Dots show mean values, bars indicate statis- 
tical dispersion of numerical data. Solid curves show power 
laws corresponding to the linear stability analysis predictions 
(see the text). The same system and simulation parameters 
as in Fig. 6. 

holds. Substituting numerical values, we get Cd = 
10 -18 Mol/cm 2 . Thus, already for very low surface con- 
centrations of proteins, effects of their diffusion can be 
neglected in the considered problem. 

Neglecting diffusion effects, the critical propulsion 
force of a single machine is estimated as 



Note that it does not depend on the protein concentra- 
tion. 

Substitution of the numerical values into equation l|18j) 
yields f c = 3 • I0~ 5 pN. Can a single molecular ma- 
chine generate the hydrodynamical propulsion force of 
that magnitude? 

Direct measurements of molecular propulsion forces 
are still not available. It is known that molecular mo- 
tors, such as myosin or kinesin, can generate mechanical 
forces of about 1 pN [T3, [H, EH , but these data refers to 
the molecules moving along microtubules and filaments, 
not to the swimmers. 

For a simple example, we consider the elementary 
three-body Purcell swimmer Its characteristic 

propulsion velocity is about V = AL/t, where AL is 
the displacement per single cycle and r is the charac- 
teristic time of the cycle (cf. [H, The viscous 
friction force of an object of linear size L that moves 
with velocity V through the liquid is by the order-of- 
magnitude fo = GirpLV. Thus, propulsion at velocity 
V = AL/t c would require the propulsion force about 
/o = 6npLAL/t c . Choosing L w 50 nm, AL = Q.1L and 
t = I ms and considering water as the liquid, we obtain 



for the molecular propulsion force the rough estimate of 
fo = I0 _3 pN. Similar estimates are obtained for other 
known elementary swimmers. Note that the actual aver- 
age propulsion force / of a machine additionally depends 
on the frequency of machine cycles and the angle with re- 
spect the interface, therefore, according to equation (JTJ) , 
this gives the estimate of the maximum propulsion force 
under the energy-supply saturation conditions and the 
most efficient orientation of the machine. 

While molecular propulsion forces are quite small, ac- 
cording the above estimates they would still be sufficient 
to induce the film instability and lead to the interface tur- 
bulence. Therefore, we conclude that the experimental 
observation of the predicted effects is principally possible. 

The experiments aimed at detecting instabilities of 
thin liquid layers induced by floating actively operating 
machines would allow to directly estimate actual propul- 
sion forces generated by particular biomoleculcs. Investi- 
gations of such noncquilibrium hydrodynamical systems 
would be very important from the perspective of active 
micro fluidics, where active motions in thin liquid layers 
are produced by floating biomolecular propellers. 
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APPENDIX A: HYDRODYNAMIC EQUATIONS 

We consider the situation when the film thickness ho 
is the smallest characteristic length of the system. In 
this case, the lubrication approximation can be used 
which corresponds to an expansion in the small parame- 
ter e = ho/X, with A being the characteristic wavelength 
of the patterns in the lateral direction. Hydrodynamical 
equations in the lubrication approximation are simplified 
and take the form (see, e.g., (l3j |) 

fid zz u = Vp, 

d z p = 0, (Al) 

with the incomprcssibility condition 

fu + d z w = 0. (A2) 

Here, w is the vertical component of the fluid velocity 
and u is its horizontal component; V is the differential 
operator acting only on the horizontal coordinates. Note 
that, as follows from these equations, pressure p is con- 
stant along the vertical direction. 

The boundary conditions at the bottom of the film, in 
contact with the solid support, are 

w = tT=0atz = 0. 
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At the free air-liquid interface z = h, the balance of hor- 
izontal and vertical forces should separately hold, imply- 
ing that 



[id z u = V<T, 

p = — a\7 2 h — p n 



(A3) 
(A4) 



where a is the surface tension coefficient and p m is the ad- 
ditional pressure produced by floating active molecules. 

Conservation of the total film volume implies moreover 
that the equation. 



numerical integration method has been employed in our 
study. Below, its brief description is provided. 

There are two nonlinear differential equations to solve, 
one for the thickness H and the other for the concentra- 
tion C of the surfactant. We define vector U = (H, C) 
and formally write both equations as 



d t U = F(U) (Bl) 
The matrix operator F can be decomposed as F = LC/ + 
N(C7) into its linear (LC/) and nonlinear (N(C/)) parts. 
They are defined as 



dth + V 




dz = 



(A5) 



should hold. 

Integrating equations. (|A1|) and taking into account the 
boundary condition (|A4j) . the horizontal velocity of the 
flow can be determined. 



hz 



zVcr 



(A6) 



The vertical velocity is then obtained by using the in- 
comprcssibility condition, 



w = -^- (\7 2 p - z 2 lij - z 2 WpWh + z 2 \J 2 a 

[ (A7) 
Substituting these expressions into equation (|A5[) and 
integrating over z, the final form of the interface equation 
is derived, 



SF(U) 



SU 



N(C/) = F(U)-LU. 



(B2) 



U=U 



To compute U(t + At) at the next time step t + At, 
the mixed semi-implicit method is employed: 



U(t + At)) - MLU{t + At) = U(t) + AtN(U{t)). (B3) 

Thus, the implicit method is used to determine the con- 
tribution from the linear part in time t + At and the 
explicit method is employed to compute the nonlinear 
part contribution. 



Applying the inverse matrix operator (1 
both sides of this equation yields 



AtL)- 1 to 



U(t + At) = (1 - AtL) U(t) + AtN(U(t)) 



Recombining some terms and using F = LC/ + N(C/), 
the final finite-difference equation is obtained 



1 - //i 3 - h 2 - \ 
^ = -v(yVp- y w). (A8) 

APPENDIX B: NUMERICAL INTEGRATION 
METHOD 

While simple explicit Euler methods are frequently em- 
ployed in numerical investigations of nonlinear reaction- 
diffusion models, such methods easily become numeri- 
cally unstable when hydrodynamical microfluidics equa- 
tions are considered. Therefore, a specially constructed 



U(t + At) = U(t) + Ai(l - AtL)" 1 F(t7(t)) 

The calculation of the inverse differential matrix op- 
erator (1 — AtL) -1 is most conveniently performed by 
transforming the equation to the Fourier space. Here, it 
is important that the term F(U(t)) is determined before 
applying the Fourier transformation. 

Our simulations using this numerical method have 
been performed only in the one-dimensional case, H = 
H{x,t) and C = C(x,t). In two dimensions, the method 
becomes complicated and such simulations have not been 
undertaken. 
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